/*********************************************************************
This is a bc script, written by Karl dahlke.
You must include the library,    bc -l sunset
Then usage is   sunset_az(latitude, month, day)
Gives the azymuth of sunset at the given latitude on the given day.
Type quit when done.

Overview:

A friend of mine in Michigan told me that sunset is, at times,
a bit north of west, not due west, but west by northwest. 🌇
The sun never gets farther than 23.5 degreees north latitude, and we are 45 degrees north latitude. How can sunset ever be north of us? Seems like a contradiction, doesn't it?

Ok I'm slowly getting a grip on this.
First, we have to remember some of our noneuclidean geometry. A "line" on the sphere is a great circle or part of a great circle. Part of a circle that is as large as it can be on the sphere.
Two lines can never be parallel, and the angles of a triangle sum to more than 180 degrees.

Hang the earth in space with the sun in front shining the front half of the earth. The back half is dark.
The line around the earth, as we look at it, is the terminator, a line that separates day and night.
It's a proper line because it is a great circle around the earth.
The earth has spun so that you, Michigan, is on the right edge. The terminator runs right through your house.
The middle of the earth is approximately Hawaii, as we have set it up.
You stand and look left, look west, and see the sunset.
Because you stand on the terminator, you are seeing the sunset, and in another ten minutes the earth turns and slides you into night, as your house goes to the back side of the earth.
So this is the setup.

Now we account for the season.
Simplest case is the equinox, the sun is directly over the equator.
As we look at the lit side of the earth, the equator runs right through the middle of it.
It is slanted up, or down, depending on which equinox, but it cuts through the middle of the earth.
The sun is a ball of light in front of the earth, floating just above the equator.
Let s be the point on the earth with the sun overhead, the point on the equator with the sun overhead.
Let p be the point where you stand, in Michigan. Remember that p is on the edge.
Let q be the point on the equator that is directly south of you.
Head due south for about 3,100 miles and there is the equator and there is q.
In our setup q is also on the edge of the earth, on the terminator.
The line from p to q runs along the terminator, like the edge of a coin.
pqs makes a triangle on the sphere.
Think of the triangle as pointing left, with s as the apex.
This is an isosceles triangle. You wouldn't think so in the plane, but it is on the sphere.
ps = qs, at about 6,200 miles.
The base of the triangle is pq, from your position south to the equator, at 3,100 miles.
It almost looks like a slice of pie.
The base angles are 90 degrees.
At q, on the equator, look straight west and go around to s.
At p, where you stand, look straight west and go around to s. Thus both right angles.
The angle at s is not zero. The three angles sum to more than 180 degrees. That's ok on the sphere.
Therefore, at equinox, either equinox, sunset is due west.
And it doesn't matter the latitude of p.
You could live south in Alabama, or north in Canada, it's still a triangle, with base from p to q, the base angles are still 90 degrees, and sunset is due west.

Now let's say it is the vernal equinox, and move into spring.
The sun appears to move northward in the sky.
It's the same picture, but the equator is down and to the right, not slicing directly through the middle.
The point s is a tad north of the equator.
It could be just above the equator, or all the way up to 23.5 degrees at solstice on June 21.
p is still where it was before, where you stand, on the terminator, on the very edge of the earth as we view it.
You are about to slip into night. You look in the general direction of west, from p to s, and there is sunset.
The equator has moved down and to the right, so q is actually behind the earth, in darkness.
Straight south from p slides behind the earth and reaches q.
This is the base of the triangle.
The triangle is not isosceles any more. qs is longer than ps.
That means you have to look just a bit above west, west by northwest, to see s, to see the sunset.
Sunset gets a little more north from March 21 to June 21, then slides back down from June 21 to September 21, when it is due west again.
The other half of the year it is south of west, west by southwest.
That's not surprising, but that sunset is actually a tad north of us for half the year, that still surprises me a bit, even after I have proven it mathematically.
This is my intuition, and I feel pretty good about it now,
so let's begin the calculations, which are at times approximations.

We need to understand s relative to the equator, on any given day of the year.
I'm going to assume earth has a circular orbit around the sun, which it doesn't, but close enough.
The north pole leans toward the sun at the start of our summer, and away from the sun in winter, and so on.
The sun is iluminating the front of the earth, as we face it.
s is always at the center of our picture.
n is the north pole.
It is at the left edge of the earth, 66.5 degrees up, at the vernal equinox.
Over the next 6 months n slides around front, and over to the right edge, always 66.5 degrees up.
The equator moves in lockstep with n.
At the start the equator runs directly through s and slants up.
When n is directly above s, the equator is below s, making something like a smile.
When n is on the right edge, the autumnal equinox, the equator passes through s again, this time slanting down and to the right.

give the earth a radius of 1.
Let the x axis point left, y points towards us through s and out to the sun, and z points up.
This is a nonstandard coordinate system, but it is convenient for this work.
Lett t = 23.5 degrees, the tilt of the earth.
n traces a circle about the z axis, at z = cos(t), with radius sin(t).
The angle w of n runs from 0 to 2π throughout the year.
The year starts at March 21, the vernal equinox, with n on the left edge.
n = sin(t)*cos(w), sin(t)*sin(w), cos(t).

From March 21 to June 21 is 92 days, then 92 days from June 21 to September 21.
Center each day at noon, time doesn't matter here.
That's 184 days which is more than half the year.
So w moves from 0 to π, and n traces half a circle, in more than half a year.
That's ok because the earth is farther from the sun in the summer, and moves a little slower in the summer, than the winter.
This approximation is a nod to the elliptical orbit,
that this half of the circle is a little slower than the other half.
So then, w is prorated by the number of days in spring and summer, or the number of days in fall and winter.
Feb 29 is not handled, sorry.
You can improve this script by adding a flag for leapyear,
then allow for feb 29.

Let me know if you want me to convert this to python or perl or other such.
*********************************************************************/

scale = 7
pi = a(1) * 4
/* switch between degrees and radians */
define d2r(u)
{
return u / 180 * pi;
}
define r2d(u)
{
return u / pi * 180;
}

/* find w from the month and day */
define find_days(m, d)
{
auto days
days = 0
if(m == 3 && d >= 21 || m > 3 && m < 9 || m == 9 && d <= 21) {
/* spring and summer */
if(m == 3) return d - 21
days = 10
if(m == 4) return days + d;
days += 30
if(m == 5) return days + d;
days += 31
if(m == 6) return days + d;
days += 30
if(m == 7) return days + d;
days += 31
if(m == 8) return days + d;
days += 31
/* at this point m is 9 */
return days + d;
}
/* fall and winder */
days = 184
if(m == 9) return days + d - 21;
days += 9;
if(m == 10) return days + d;
days += 31;
if(m == 11) return days + d;
days += 30;
if(m == 12) return days + d;
days += 31;
if(m == 1) return days + d;
days += 31;
if(m == 2) return days + d;
days += 28;
/* at this point m is 3 */
return days + d;
}

define find_w(m, d)
{
auto days
days = find_days(m, d)
if(days <= 184) return days / 184 * pi
return (days-184) / 181 * pi + pi
}

/*********************************************************************
Start at n, head 90 degrees in any direction,
turn 90 degrees, and run along the equator.
We may as well run from n through s.
s is at 0,1,0
The cosine of the angle between any two unit vectors is their dot product.
s.n = sin(t) * sin(w)
the arc cosine of this is the angle from n to s.
Subtract from 90 degrees and get the distance from s to the equator.
arc cosine is not native to bc, we have to write it.
If you switch to another language, and call acos directly,
be careful, it might not do what you want on a negative argument.
acos(-0.01) should equal approx π/2 + 0.01.
*********************************************************************/

define acos(u)
{
if(u == 0) return pi/2
if(u > 0) return a(sqrt(1-u^2)/u)
/* u is negative, it's fall or winter */
return pi - a(sqrt(1-u^2)/(-u))
}

/* tilt of the earth. Later on we will use t for tangent function,
but don't worry, in bc the name space for variables
is separate from the name space for functions.
This is not typical of other languages. */
t = d2r(23.5)

define sun_above(w)
{
return pi/2 - acos(s(t)*s(w))
}

/*********************************************************************
Now shift the picture.
s is not in the center any more.
The equator runs across the center of the earth, left to right,
and cuts the earth in half. The sun
shines down at a point called s, which floats above or below the equator.
p is the point where you stand, on the right edge of the earth.
Designate p by its latitude, that's all you need.
Latitudes at or above the arctic circle give nonsensical results.
The sun might not set at all; who knows.
Best to stay below 60 degrees.
We know the latitude of s, from the sun_above function.
We need the longitude of s.
That's what we compute next, slong, the longitude of s.
I'm measuring longitude from left to right, 0 to 180.
This is x = 1 to -1 in our coordinate system.
It is sunset at p, so s to p is 90 degrees.
The dot product of the two vectors of s and p is 0.
plat is the latitude of p in radians.
The p vector is -cos(plat), 0, sin(plat).
the s vector is cos(slong)*cos(above), sin(slong)*cos(above), sin(above).
cos(plat) * cos(slong) * cos(above) = sin(plat) * sin(above)
cos(slong) = tan(plat) * tan(above)
slong = acos(tan(plat) * tan(above))
First define tangent, which is not native to bc.
*********************************************************************/

define t(u) { return s(u) / c(u) }

define find_slong(plat, above)
{
return acos(t(plat) * t(above))
}

/*********************************************************************
There is a plane determined by p, s, and the center of the earth.
When you stand at p and look toward the sunset you are confined to this plane.
Let's find the equation of this plane.
Because the plane contains the origin, it has the form ax + by + cz = 0.
restrict to y = 0, and determine c = 1 and a = tan(plat).
Go back to the s vector, which was a bit ugly.
cos(slong)*cos(above), sin(slong)*cos(above), sin(above).
a * cos(slong)*cos(above) + b * sin(slong)*cos(above) + c * sin(above) = 0
tan(plat) * cos(slong)*cos(above) + b * sin(slong)*cos(above) + sin(above) = 0
sin(slong) and cos(above) can never be 0.
tan(plat) * cos(slong) + b * sin(slong) + tan(above) = 0
b = - (tan(plat) * cos(slong) + tan(above)) / sin(slong)
Now we have our plane ax + by + cz = 0.
Can this plane ever contain the y axis?
for any p, s sits on the equator at the equinox.
This sets slong = 90 degrees.
On any other day, s is off the equator, and slong changes to adjust.
Our plane contains the y axis only on the equinox.
In cofirmation, b = 0 only for the equinox.
We will write some exception code: if equinox, look due west,
otherwise, run some calculations, where we know b is nonzero.
You want to stand at p and look toward s.
You want to look within the plane that contains s p and origin.
You also have to look within the plane that is tangent to the earth, the horizon.
The intersection of two planes is a line, the look vector, which I will call l.
l is constrained by a b and c.
The vector sin(plat),0,cos(plat) points due north from p.
This is a unit vector. Call this vector n.
The vector -cos(plat),0,sin(plat) points straight up from the ground.
This is a unit vector. Call this vector up.
l.up = 0.
we have two equations which will define l.
a*lx + b*ly + c*lz = 0
-cos(plat)*lx + sin(plat)*lz = 0
lx = lz * tan(plat).
a*lz*tan(plat) + b*ly + c*lz = 0
lz*tan(plat)^2 + b*ly + lz = 0
lz*(tan(plat)^2 + 1) + b*ly = 0
ly = - lz*(tan(plat)^2 + 1) / b.
Set lz = 1 and compute lx and ly.
We have our look vector.
The azymuth is now the angle between l and n.
this is the arc cosine of l.n divided by the lengths of the two vectors.
n is already a unit vector so don't worry about that.
acos(l.n / |l|)
Here we go.
*********************************************************************/

define sunset_az(p_degrees, m, d)
{
auto plat, w, above, slong, a, b, lx, ly, lz, ll, az
/* turn p latitude degrees into radians */
plat = d2r(p_degrees)
/* assume month and day are valid */
w = find_w(m, d)
above = sun_above(w)
/* if equinox, the answer is always due west */
if(above == 0) return 270
/* another special case, standing on the equator */
if(plat == 0) return 270 + r2d(above);
/* more calculations needed. Shift picture so equator is the center.
Latitude of s is above, find its longitude. */
slong = find_slong(plat, above)
/* coefficients for the plane, c is always 1 so we don't need that */
a = t(plat)
b = - (t(plat) * c(slong) + t(above)) / s(slong)
/* look vector */
lz = 1
/* I'm not going to include lz over and over in the formulas */
lx = t(plat)
ly = - (lx^2 + 1) / b
ll = sqrt(lx^2 + ly^2 + 1)
/* l dot n */
az = s(plat)*lx + c(plat)
/* make sure l wasn't pointing east */
if(ly < 0) az = -az;
/* find the angle */
az = acos(az / ll)
az = 2*pi - az
az = r2d(az)
return az
}